library(tidyverse)
library(lme4)

options(scipen=999)
df = read.csv(here::here('avg_tensor_by_roi_wide.csv'),
              colClasses = c('subject' = 'factor',
                             'site' = 'factor',
                             'visit' = 'factor'))
outcomes = df %>% 
  select(where(is.numeric)) %>% 
  colnames

Summary

The purpose of the DTIMSCAMRAS study is to test for site effects on DTI metrics of different brain regions in the context of a traveling subjects design with a harmonized imaging protocol. Participants (n=11) with stable Multiple Sclerosis (MS) were scanned in 4 different locations: National Institute of Health (NIH), Brigham and Women’s Hospital (BWH), Johns Hopkins University (Hopkins), and The University of Pennsylvania (Penn). Two scans were collected per site visit, and participants were re-positioned between scans. Penn, BWH and NIH used Siemens scanners, while Hopkins used a Phillips scanner. The imaging modalities that were collected include T1s (MPRAGE), FLAIR, and Diffusion-weighted images (DWI).

Preprocessing

To preprocess imaging data, qsiprep was used for bias correction, skull-stripping and co-registration of the T1s and DWIs, and specialized denoising and artifact correction for the DWIs. Then, DTIFIT was used to compute the following scalar maps:

  • Fractional Anisotropy
  • Mean Diffusivity
  • Axial Diffusivity
  • Radial Diffusivity

Concurrently, several segmentation pipelines were used on the T1s (+ FLAIRs in the case of MIMOSA) with the purpose of defining particular regions of interest (ROIs):

  • White Matter and Gray Matter Segmentation methods:
    • ANTs ATROPOS
    • FSL FAST
    • Multi-Atlas Segmentation with Joint Label Fusion using the MUSE templates (JLF WM and JLF GM)
  • Thalamus Segmentation
    • FSL FIRST
    • Multi-Atlas Segmentation with Joint Label Fusion using the OASIS templates (JLF THALAMUS)
  • Lesion Segmentation methods:
    • MIMOSA

To remove lesions appearing in GM or WM (or thalamus) from the calculation, mimosa masks were inverted and applied to the segmentation label maps in order to zero out the tissue masks at the location of each lesion. These adjusted segmentation results were then used to average the DTI scalar maps across all voxels belonging to each of the ROIs. The 36 resulting DTI metrics make up the set of outcome variables which are the focus of the site effects test. Below, the outcome variables are listed.

data.frame(OUTCOME = outcomes) %>% 
  separate(OUTCOME, c("SEGMENTATION", "ROI", "SCALAR"), remove = FALSE)

Note the general format: for example, the ATROPOS_GM_AD variable denotes the average Axial Diffusivity (AD) across the Gray Matter (GM) regions defined by the Atropos Segmentation.

Testing site effects

Testing for site effects was performed using two approaches: a permutation-based approach and a deviance test of mixed models.

Additionally, as Hopkins data was acquired using a Philips scanner, which is in contrast to all other sites, the two testing approaches were performed again while excluding Hopkins data.

Permutation-Based F-test

A permutation-based F-test was used to test for site effects. For each subject, the squared difference between all possible pairs of intra-site and inter-site measures were computed. For instance, take the first row of the subject data.frame below, which has BWH as the “reference site”: the squared differences between that row and the 6 non-BWH rows are computed for \(y\)——this process is repeated for all rows (and the results averaged) to arrive to the mean inter-site squared difference for this subject; to calculate the mean intra-site squared differences, the squared differences of the 4 intrasite pairs are averaged.

df %>% 
  filter(subject == '01001') %>% 
  select(!where(is.numeric), ATROPOS_GM_AD) %>% 
  rename(y = ATROPOS_GM_AD)

These differences were averaged across all subjects resulting in the Mean Squared difference Between Sites (\(\overline{MS_{B}}\)) and the Mean Squared difference Within Sites (\(\overline{MS_{W}}\)). The test statistic is composed of their ratio:

\[ F = \frac{\overline{MS_{B}}}{\overline{MS_{W}}} \]

For each metric, the observed test statistic was compared to a distribution of test statistics derived from 10000 permutations (i.e. a null distribution). The site labels were permuted within each subject. The proportion of null results higher than the observed test statistic served as a preliminary p-value \(p_0\). To prevent p-values of 0 (when the observed statistic was higher than all permuted results), the p-value was shifted using the following formula:

\[ p = \frac{10000*p_0 + 1}{10000 + 1}\]

Crossed-design Mixed Linear Models

As an additional test of site effects, mixed models were performed corresponding to a crossed design using the lmer function of the lme4 package. For each outcome variable, two mixed-effects models were fitted: The formula for the base model includes a random intercept for site as well as a random intercept for subject. The formula for the extended model also includes these terms in addition to a fixed effect term for site.

Base Model:

y ~ (1|subject) + (1|site)

Extended model:

y ~ site + (1|subject) + (1|site)

For each outcome, these two models were compared through a deviance test using the anova function in R. P-values for these tests are reported as a test of site effects.

Source code for this project can be found on GitHub.

Key Findings

  • Permutation tests revealed significant site effects in 31 out of 36 outcomes after FDR adjustment when all data is included. After excluding Hopkins data, the number of significant site effects decreased to 13. Most of the DTI metrics (i.e. outcomes) that showed differences were derived from thalamus segmentation results.
  • For mixed models, most outcome variables (13 out of 36) show a significant test after FDR adjustment indicating the presence of site effects when all data is included. When Hopkins data is excluded, none of the outcome variables show site effects.

Questions for follow up

  • It is unclear to what extent potential segmentation differences across sites might contribute to site effects, i.e. to what extent do site differences stem from ROI coverage differences (due to scanner effects on the T1s), and to what extent are the differences due to scanner effects on DWI images?

Data Inventory

Subjects 03001 and 03002 did not complete imaging at all 4 sites. Subject 03001 completed imaging at Penn and BWH only and 03002 completed imaging at BWH, Hopkins, and the NIH only.

Subject 02001 is missing DTI data from their NIH visit.

Scans at site per subject

t_df <- df %>% 
   count(subject, site, .drop = FALSE) %>% 
  tidyr::pivot_wider(names_from = site, values_from = n) %>% 
  mutate(`Row Totals` = BWH + Hopkins + NIH + Penn) 
t_df <- t_df %>% 
  summarize(across(where(is.integer), sum)) %>% 
  bind_rows(t_df, .) %>% 
  mutate(across(where(is.factor), as.character))
t_df[is.na(t_df)] <- "Column Totals"
t_df

For Atropos (substracting failed segmentations)

Atropos segmentation failed for 6 images: 04001NIH01, 01003NIH01, 03002NIH02, 04003NIH01, 04001BWH02, and 03001BWH01.

knitr::include_graphics(here::here(c('misc/atropos_success.png', 'misc/atropos_fail.png')))
Left: Example of successful Atropos segmentation; Right: Example of a failed Atropos segmentation with one ROI missingLeft: Example of successful Atropos segmentation; Right: Example of a failed Atropos segmentation with one ROI missing

Left: Example of successful Atropos segmentation; Right: Example of a failed Atropos segmentation with one ROI missing

The following table excludes the failed atropos segmentations. Note subject 03001 only has 3 images after excluding failed atropos segmentations.

t_df <- df %>% 
  unite(sub, subject, site, visit, sep = '-') %>%
  filter(!sub %in% c('04001-NIH-01', '01003-NIH-01', '03002-NIH-02', '04003-NIH-01', '04001-BWH-02', '03001-BWH-01')) %>%
  separate(sub, c('subject', 'site', 'visit')) %>%
  mutate_if(is.character, as.factor) %>% 
  count(subject, site, .drop = FALSE) %>% 
  tidyr::pivot_wider(names_from = site, values_from = n) %>% 
  mutate(`Row Totals` = BWH + Hopkins + NIH + Penn)

t_df <- t_df %>% 
  summarize(across(where(is.integer), sum)) %>% 
  bind_rows(t_df, .) %>% 
  mutate(across(where(is.factor), as.character))
t_df[is.na(t_df)] <- "Column Totals"
t_df

The failed atropos segmentations are included in the following visualizations but were excluded for the site effects tests.

TODO: compare harmonization with and without hopkins

Visualizations

Raw Data

Fractional Anisotropy

Densities

Densities are colored by site, while the black line represents the overall density aggregated across sites. Distribution of the different sites tend to cluster together except for JLF Thalamus. Additionally, note the long tails caused by outliers in the ATROPOS distributions.

plot_densities <- function(var){
  df %>% 
    ggplot(aes_string(x=var, color='site')) + 
    geom_density() + 
    stat_density(aes_string(x = var), geom = "line", color = "black") +
    theme_bw() + 
    xlab(str_replace_all(var, "_", " "))
}

vars <-df %>% 
  select(ends_with('FA')) %>% 
  colnames()
dplots <- lapply(vars, plot_densities) %>% 
  setNames(vars)
for (name in names(dplots)){
  cat("#####",  str_replace_all(name, "_", " "), "\n")
  print(dplots[[name]])
  cat("\n\n")
}
ATROPOS GM FA

ATROPOS WM FA

FAST GM FA

FAST WM FA

FIRST THALAMUS FA

JLF GM FA

JLF WM FA

JLF THALAMUS FA

MIMOSA LESION FA

Box plots

Box plots show difference in medians across sites for most variables.
Additionally, note the severe outliers in the ATROPOS metrics which correspond to the failed segmentations. For FIRST THALAMUS, subject 01002 could be considered an outlier; their FA values are consistent across sites and scans, however, suggesting this variation is meaningful and that this subject should be retained.

plot_avg_tensor_values <- function(tensor){
  df %>% 
    select(!is.numeric, ends_with(tensor)) %>% 
    pivot_longer(ends_with(tensor), names_to = 'seg', values_to = 'average') %>% 
    separate(seg, into = c('segmentation', 'roi', 'tensor')) %>% 
    unite(segmentation, segmentation, roi, sep = " ") %>% 
    mutate(segmentation = factor(segmentation,
                                 levels = c("ATROPOS GM", "ATROPOS WM", "FAST GM", "FAST WM",
                                            "JLF GM", "JLF WM", "JLF THALAMUS",  "FIRST THALAMUS", "MIMOSA LESION"))) %>%
    ggplot(aes(x=site, y=average)) + 
    geom_boxplot(outlier.shape = NA) + 
    geom_jitter(aes(color = subject, shape=visit), alpha=0.5) +
    #facet_grid(site~.) + 
    facet_grid(segmentation ~ .) +
    coord_flip() +
    ggtitle(sprintf("Mean %s in ROI by site", tensor)) +
    theme_bw() + 
    # theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust=1)) +
    xlab("") + 
    ylab("") + 
    theme(strip.text.y = element_text(angle = 0)) + 
    scale_x_discrete(expand = c(0.15, 0.15))
}

plot_avg_tensor_values('FA')

Mean Diffusivity

Densities

For Mean Diffusivity, distributions vary more widely across sites. In particular, MIMOSA LESION MD values are quite different for Hopkins data.

vars <-df %>% 
  select(ends_with('MD')) %>% 
  colnames()
dplots <- lapply(vars, plot_densities) %>% 
  setNames(vars)
for (name in names(dplots)){
  cat("#####", str_replace_all(name, "_", " "), "\n")
  print(dplots[[name]])
  cat("\n\n")
}
ATROPOS GM MD

ATROPOS WM MD

FAST GM MD

FAST WM MD

FIRST THALAMUS MD

JLF GM MD

JLF WM MD

JLF THALAMUS MD

MIMOSA LESION MD

Boxplot

As before, medians appear to be different across sites.

plot_avg_tensor_values('MD')

Radial Diffusivity

Densities

Radial Diffusivity also shows some variations across sites (perhaps less than MD). JLF GM RD in NIH shows multiple peaks, and JLF WM RD shows skewed distributions for all sites. As before the metric for MIMOSA shows a different distribution for Hopkins data.

vars <-df %>% 
  select(ends_with('RD')) %>% 
  colnames()
dplots <- lapply(vars, plot_densities) %>% 
  setNames(vars)
for (name in names(dplots)){
  cat("#####",  str_replace_all(name, "_", " "), "\n")
  print(dplots[[name]])
  cat("\n\n")
}
ATROPOS GM RD

ATROPOS WM RD

FAST GM RD

FAST WM RD

FIRST THALAMUS RD

JLF GM RD

JLF WM RD

JLF THALAMUS RD

MIMOSA LESION RD

Boxplot

As before, medians appear to be differet across sites.

plot_avg_tensor_values('RD')

Axial Diffusivity

Densities

For Axial Diffusivity, Hopkins distributions show marked differences across most ROIs compared to other sites. The difference is quite pronounced for MIMOSA.

vars <-df %>% 
  select(ends_with('AD')) %>% 
  colnames()
dplots <- lapply(vars, plot_densities) %>% 
  setNames(vars)
for (name in names(dplots)){
  cat("#####",  str_replace_all(name, "_", " "), "\n")
  print(dplots[[name]])
  cat("\n\n")
}
ATROPOS GM AD

ATROPOS WM AD

FAST GM AD

FAST WM AD

FIRST THALAMUS AD

JLF GM AD

JLF WM AD

JLF THALAMUS AD

MIMOSA LESION AD

Boxplot

As before, medians show differences, particularly for Hopkins.

plot_avg_tensor_values('AD')

Harmonized Data

Combat estimates

# load harmonized data
harmonized_df <- read.csv(here::here('results/avg_tensor_by_roi_wide_harmonized.csv'),
                          colClasses = c('subject' = 'character',
                             'site' = 'character',
                             'visit' = 'character'))
harmonized_no_hopkins_df <- read.csv(here::here('results/avg_tensor_by_roi_wide_no_Hopkins_harmonized.csv'))

# load models
combat_model <- readRDS(here::here('results/avg_tensor_by_roi_wide_model.csv'))[c("gammahat", "delta2hat", "gammastarhat", "delta2starhat")]
combat_model_no_hopkins <- readRDS(here::here('results/avg_tensor_by_roi_wide_no_Hopkins_model.csv'))[c("gammahat", "delta2hat", "gammastarhat", "delta2starhat")]

# model parameters
# what's the best way to plot combat estimates

Fractional Anisotropy

Densities

Densities are colored by site, while the black line represents the overall density aggregated across sites. Distribution of the different sites tend to cluster together except for JLF Thalamus. Additionally, note the long tails caused by outliers in the ATROPOS distributions.

plot_densities <- function(var){
  harmonized_df %>% 
    ggplot(aes_string(x=var, color='site')) + 
    geom_density() + 
    stat_density(aes_string(x = var), geom = "line", color = "black") +
    theme_bw() + 
    xlab(str_replace_all(var, "_", " "))
}

vars <- harmonized_df %>% 
  select(ends_with('FA')) %>% 
  colnames()
dplots <- lapply(vars, plot_densities) %>% 
  setNames(vars)
for (name in names(dplots)){
  cat("#####",  str_replace_all(name, "_", " "), "\n")
  print(dplots[[name]])
  cat("\n\n")
}
FAST GM FA

FAST WM FA

FIRST THALAMUS FA

JLF GM FA

JLF WM FA

JLF THALAMUS FA

MIMOSA LESION FA

ATROPOS GM FA

ATROPOS WM FA

Box plots

Box plots show difference in medians across sites for most variables.
Additionally, note the severe outliers in the ATROPOS metrics which correspond to the failed segmentations. For FIRST THALAMUS, subject 01002 could be considered an outlier; their FA values are consistent across sites and scans, however, suggesting this variation is meaningful and that this subject should be retained.

plot_avg_tensor_values <- function(tensor){
  harmonized_df %>% 
    select(!is.numeric, ends_with(tensor)) %>% 
    pivot_longer(ends_with(tensor), names_to = 'seg', values_to = 'average') %>% 
    separate(seg, into = c('segmentation', 'roi', 'tensor')) %>% 
    unite(segmentation, segmentation, roi, sep = " ") %>% 
    mutate(segmentation = factor(segmentation,
                                 levels = c("ATROPOS GM", "ATROPOS WM", "FAST GM", "FAST WM",
                                            "JLF GM", "JLF WM", "JLF THALAMUS",  "FIRST THALAMUS", "MIMOSA LESION"))) %>%
    ggplot(aes(x=site, y=average)) + 
    geom_boxplot(outlier.shape = NA) + 
    geom_jitter(aes(color = subject, shape=visit), alpha=0.5) +
    #facet_grid(site~.) + 
    facet_grid(segmentation ~ .) +
    coord_flip() +
    ggtitle(sprintf("Mean %s in ROI by site", tensor)) +
    theme_bw() + 
    # theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust=1)) +
    xlab("") + 
    ylab("") + 
    theme(strip.text.y = element_text(angle = 0)) + 
    scale_x_discrete(expand = c(0.15, 0.15))
}

plot_avg_tensor_values('FA')

Mean Diffusivity

Densities

For Mean Diffusivity, distributions vary more widely across sites. In particular, MIMOSA LESION MD values are quite different for Hopkins data.

vars <- harmonized_df %>% 
  select(ends_with('MD')) %>% 
  colnames()
dplots <- lapply(vars, plot_densities) %>% 
  setNames(vars)
for (name in names(dplots)){
  cat("#####", str_replace_all(name, "_", " "), "\n")
  print(dplots[[name]])
  cat("\n\n")
}
FAST GM MD

FAST WM MD

FIRST THALAMUS MD

JLF GM MD

JLF WM MD

JLF THALAMUS MD

MIMOSA LESION MD

ATROPOS GM MD

ATROPOS WM MD

Boxplot

As before, medians appear to be different across sites.

plot_avg_tensor_values('MD')

Radial Diffusivity

Densities

Radial Diffusivity also shows some variations across sites (perhaps less than MD). JLF GM RD in NIH shows multiple peaks, and JLF WM RD shows skewed distributions for all sites. As before the metric for MIMOSA shows a different distribution for Hopkins data.

vars <-harmonized_df %>% 
  select(ends_with('RD')) %>% 
  colnames()
dplots <- lapply(vars, plot_densities) %>% 
  setNames(vars)
for (name in names(dplots)){
  cat("#####",  str_replace_all(name, "_", " "), "\n")
  print(dplots[[name]])
  cat("\n\n")
}
FAST GM RD

FAST WM RD

FIRST THALAMUS RD

JLF GM RD

JLF WM RD

JLF THALAMUS RD

MIMOSA LESION RD

ATROPOS GM RD

ATROPOS WM RD

Boxplot

As before, medians appear to be differet across sites.

plot_avg_tensor_values('RD')

Axial Diffusivity

Densities

For Axial Diffusivity, Hopkins distributions show marked differences across most ROIs compared to other sites. The difference is quite pronounced for MIMOSA.

vars <-harmonized_df %>% 
  select(ends_with('AD')) %>% 
  colnames()
dplots <- lapply(vars, plot_densities) %>% 
  setNames(vars)
for (name in names(dplots)){
  cat("#####",  str_replace_all(name, "_", " "), "\n")
  print(dplots[[name]])
  cat("\n\n")
}
FAST GM AD

FAST WM AD

FIRST THALAMUS AD

JLF GM AD

JLF WM AD

JLF THALAMUS AD

MIMOSA LESION AD

ATROPOS GM AD

ATROPOS WM AD

Boxplot

As before, medians show differences, particularly for Hopkins.

plot_avg_tensor_values('AD')

Permutation Testing of Site Effects

# replace ATROPOS measures with NA for select images (segmentation failed)
fill_na_atropos <- function(df){
  atropos_cols <- df %>% 
  select(contains('ATROPOS')) %>% 
  colnames()
df[(df$subject == '04001' & df$site == 'NIH' & df$visit == '01') |
   (df$subject == '01003' & df$site == 'NIH' & df$visit == '01') |
   (df$subject == '03002' & df$site == 'NIH' & df$visit == '02') |
   (df$subject == '04003' & df$site == 'NIH' & df$visit == '01') |
   (df$subject == '04001' & df$site == 'BWH' & df$visit == '02') |
   (df$subject == '03001' & df$site == 'BWH' & df$visit == '01'), atropos_cols] <- NA
  return(df)
}

df <- fill_na_atropos(df)
harmonized_df <- fill_na_atropos(harmonized_df)

The permutation-based F-test was carried out as previously described. The observed F-statistic for each of the metrics is plotted below as a black dot, while the distribution of permuted test statistics is shown as a white violin plot. All metrics showed significant site effects.

ratio_stats = readRDS(here::here("results/avg_tensor_by_roi_wide_F_df.rds"))
null_dist = readRDS(here::here("results/avg_tensor_by_roi_wide_null_F_df.rds")) %>% 
  mutate(outcome = map_chr(outcome, unique)) %>% 
  select(-rowname)

All sites (raw data)

F Distributions

Plot
plot_F_dist <- function(ratio_stats, null_dist){
  null_dist %>% 
  select(outcome, F) %>% 
  unnest() %>% 
  ggplot(aes(x=outcome, y=F)) + 
  geom_violin(draw_quantiles = c(0.95)) + 
  coord_flip() + 
  geom_point(aes(x=outcome, y=F), 
             data=tibble::rownames_to_column(ratio_stats)) 
}

plot_F_dist(ratio_stats, null_dist)

Table
ratio_stats %>% 
  arrange(desc(F))

p-values

#' get p value of F by comparing to permutation distribution
test_ratio_stat = function(ratio.stats, null.dists){
  
  p.vals = c()
  for (col in colnames(ratio.stats)){
    ratio.stat = pull(ratio.stats, col)
    null.dist = pull(null.dists, col)
    
    percentile = ecdf(null.dist)
    p.val <- (10000* (1 - percentile(ratio.stat)) + 1)/(1+10000)
    p.vals = c(p.vals, p.val)
  }
  names(p.vals) = colnames(ratio.stats)
  return(p.vals)
}

#' make a data.frame with p values
make_perm_df <- function(ratio_stats, null_dist){
  observed_F <- ratio_stats  %>%
  select(outcome, F) %>% 
  pivot_wider(names_from = 'outcome', values_from = 'F')


  null_Fs <- null_dist %>% 
    select(outcome, F) %>% 
    unnest() %>% 
    pivot_wider(names_from = 'outcome', values_from = 'F') %>% 
    unnest() %>% 
    select(!!!colnames(observed_F)) # re-order columns to match observed
  
  p_df <- data.frame(p_val = test_ratio_stat(observed_F, null_Fs))
  
  perm_df <- cbind(
    ratio_stats %>% 
      select(F),
    p_df %>% 
      mutate(pval_fdr = p.adjust(p_val, method = 'fdr'))
  )
  return(perm_df)
}

perm_df <- make_perm_df(ratio_stats, null_dist)
Plot
plot_p_perm <- function(perm_df){
  perm_df %>% 
    tibble::rownames_to_column('outcome') %>% 
    mutate(outcome = fct_reorder(outcome, pval_fdr, .desc = TRUE)) %>% 
    ggplot(aes(x=outcome, y=pval_fdr)) + 
    geom_point(size = 2) +
    geom_hline(yintercept = 0.05, color = 'red') + 
    coord_flip() + 
    ggtitle("P-values of Permutation Test") +
    ylim(0,1)
}

plot_p_perm(perm_df)

Table
perm_df %>% 
  arrange(desc(p_val))

There were 35 out of 36 significant site effects according to permutation test.

No Hopkins (raw)

Image acquisition at Hopkins was performed with a Philips scanner and distortion correction during on-scanner reconstruction. As the scanner manufacturer and reconstructions method differed than other sites, here we report the permutation results based on the non-Hopkins data only.

ratio_stats = readRDS(here::here("results/avg_tensor_by_roi_wide_no_Hopkins_F_df.rds"))
null_dist = readRDS(here::here("results/avg_tensor_by_roi_wide_no_Hopkins_null_F_df.rds"))  %>% 
  mutate(outcome = map_chr(outcome, unique)) %>% 
  select(-rowname)

F Distributions

Plot
plot_F_dist(ratio_stats, null_dist)

Table
ratio_stats %>% 
  arrange(desc(F))

p-values

perm_df <- make_perm_df(ratio_stats, null_dist)
Plot
plot_p_perm(perm_df)

Table
perm_df %>% 
  arrange(desc(p_val))

All site (harmonized)

ratio_stats = readRDS(here::here("results/avg_tensor_by_roi_wide_harmonized_F_df.rds"))
null_dist = readRDS(here::here("results/avg_tensor_by_roi_wide_harmonized_null_F_df.rds")) %>% 
  mutate(outcome = map_chr(outcome, unique)) %>% 
  select(-rowname)

F Distributions

Plot
plot_F_dist(ratio_stats, null_dist)

Table
ratio_stats %>% 
  arrange(desc(F))

p-values

perm_df <- make_perm_df(ratio_stats, null_dist)
Plot
plot_p_perm(perm_df)

Table
perm_df %>% 
  arrange(desc(p_val))

No Hopkins (harmonized)

ratio_stats = readRDS(here::here("results/avg_tensor_by_roi_wide_no_Hopkins_harmonized_F_df.rds"))
null_dist = readRDS(here::here("results/avg_tensor_by_roi_wide_no_Hopkins_harmonized_null_F_df.rds")) %>% 
  mutate(outcome = map_chr(outcome, unique)) %>% 
  select(-rowname)

F Distributions

Plot
plot_F_dist(ratio_stats, null_dist)

Table
ratio_stats %>% 
  arrange(desc(F))

p-values

perm_df <- make_perm_df(ratio_stats, null_dist)
Plot
plot_p_perm(perm_df)

Table
perm_df %>% 
  arrange(desc(p_val))

Mixed Models

All Sites (raw data)

For each of the 36 outcome variables, two mixed models (a base and an extended model) were fit and the extended model of interest was compared against the null/base model using a likelihood ratio test. The resulting p-values are reported. ICCs based on the extented model are also reported.

Of note, two of the models estimated a 0 variance for the site random effect. P-values and ICCs are not provided for these two models.

# Full models
run_full_models <- function(data){
  models <- outcomes %>% 
   lapply(function(x) lmer(reformulate("(1|subject) + (1|site)", x), data=data)) %>% 
   setNames(outcomes)
}

# Null models
run_null_models <- function(data){
  models_0 <- outcomes %>% 
    lapply(function(x) lmer(reformulate("(1|subject)", x), data=data)) %>% 
    setNames(outcomes)
}

#' Test variance term
#' @param mod full model
#' @param mod_0 null model
test_variance_term <- function(mod, mod_0){
  
  A = lme4:::anova.merMod(mod, mod_0)
  lrt_stat = as.numeric(-2*(logLik(mod_0) - logLik(mod)))
  
  #if(near(lrt_stat, A$Chisq[2])) return(NA)
  
  # Naive p-value given the chi_sq(df = 1) null distribution
  pval_naive = pchisq(q = lrt_stat, df = 1, lower.tail = FALSE) 
  
  # p-value given a 50:50 mixture of chisq(1) and chisq(2),
  # according to Self & Liang 1987
  pval_SL = pchisq(q = lrt_stat, df = 1, lower.tail = FALSE)/2 +
    pchisq(q = lrt_stat, df = 2, lower.tail = FALSE)/2
  
  return(data.frame(pval_naive = pval_naive,
                    pval_SL = pval_SL))
}

# Which random effect variance is 0 
get_zero_var_model_names <- function(models){
  get_var <- function(model) summary(model)$varcor$site[1,1]
  
  zero_var_models <- models %>% 
    purrr::map_dbl(get_var) %>% 
    as.data.frame() %>% 
    setNames('variance') %>% 
    rownames_to_column(var = 'outcome') %>% 
    filter(near(variance, 0, .Machine$double.eps)) %>% 
    pull(outcome)
  return(zero_var_models)
}

models <- run_full_models(df)
models_0 <- run_null_models(df)
zero_var_model_names <- get_zero_var_model_names(models)
See models with 0 variance
for (name in zero_var_model_names){
  cat('###', name, '\n\n')
  print(summary(models[[name]]))
  cat("\n\n")
}
## ### MIMOSA_LESION_FA 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: MIMOSA_LESION_FA ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -353.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5395 -0.3266  0.0607  0.3375  4.0443 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  subject  (Intercept) 0.0013879 0.03725 
##  site     (Intercept) 0.0000000 0.00000 
##  Residual             0.0004207 0.02051 
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.32571    0.01147   28.39
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

p-values

make_SL_df <- function(models, models_0){
  purrr::map2(models, models_0, test_variance_term) %>% 
    setNames(outcomes) %>% 
    bind_rows(.id = 'outcome') %>% 
    #filter(!(outcome %in% zero_var_model_names)) %>% 
    select(-pval_naive) %>% 
    rename(pval = pval_SL) %>% 
    mutate(pval_fdr = p.adjust(pval, method = 'fdr')) %>% 
    arrange(pval)
}

p_df <- make_SL_df(models, models_0)
Plot
plot_p_SL <- function(p_df){
  p_df %>% 
    mutate(outcome = fct_reorder(outcome, pval_fdr, .desc = TRUE)) %>% 
    ggplot(aes(x=outcome, y=pval_fdr)) + 
    geom_point(size = 2) +
    geom_hline(yintercept = 0.05, color = 'red') + 
    coord_flip() + 
    ggtitle("P-values of Likelihood Ratio Test (Self & Liang 1987)") +
    ylim(0,1)
}

plot_p_SL(p_df)

Table
options(scipen = 999)
p_df

ICC

ICC for the 35 remaining models is shown

options(scipen=999)
#' Get ICC for a single model
get_icc <- function(df){

  outcomes %>% 
    map(~lmer(reformulate("(1|subject) + (1|site)", .x), data=df)) %>% 
    purrr::map(~ performance::icc(.x, by_group = TRUE) %>% as.data.frame) %>% 
    setNames(outcomes) %>% 
    bind_rows(.id = "outcome") %>% 
    select_if(~ !all(is.na(.))) %>% # eliminate column of NAs
    na.omit() %>% # eliminate rows with NA
    pivot_wider(names_from = 'Group', values_from = 'ICC', names_prefix = 'ICC_') %>% 
    relocate(ICC_site, .after = 'outcome') %>%
    arrange(desc(ICC_site))
}
# extract observed icc_df
(icc_df <- get_icc(df))
Site
## Make violin plots
## AND make line
plot_ICC_site <- function(icc_df, boot_df){
  icc_df %>% 
    left_join(boot_df %>% 
                select(-ICC_site, -ICC_subject), by = 'outcome') %>% 
    mutate(outcome = fct_reorder(outcome, ICC_site)) %>% 
    ggplot(aes(x=outcome, y=ICC_site)) + 
    geom_errorbar(aes(ymin=lwr_site, ymax=uppr_site)) +
    geom_point() + 
    coord_flip() + 
    ggtitle("ICC of Site Random Effect") + 
    ylim(0,1)
}

boot_df <- readRDS(here::here('results/avg_tensor_by_roi_wide_boot_icc.rds'))
plot_ICC_site(icc_df, boot_df)

Subject
plot_ICC_subject <- function(icc_df, boot_df){
  icc_df %>% 
    left_join(boot_df %>% 
              select(-ICC_site, -ICC_subject), by = 'outcome') %>% 
    mutate(outcome = fct_reorder(outcome, ICC_subject)) %>% 
    ggplot(aes(x=outcome, y=ICC_subject)) + 
    geom_errorbar(aes(ymin=lwr_subject, ymax=uppr_subject)) +
    geom_point() + 
    coord_flip() + 
    ggtitle("ICC of Subject Random Effect") + 
    ylim(0,1)
}
plot_ICC_subject(icc_df, boot_df)

After FDR-adjustment, 34 of the 36 variables show site effects according to the mixed models.

No Hopkins (raw data)

Here, we report results from the mixed models carried out only on the non-Hopkins data.

models <- df %>% 
  filter(site != 'Hopkins') %>% 
  run_full_models()

models_0 <- df %>% 
  filter(site != 'Hopkins') %>% 
  run_null_models()

zero_var_model_names <- get_zero_var_model_names(models)
See models with 0 variance
for (name in zero_var_model_names){
  cat('###', name, '\n\n')
  print(summary(models[[name]]))
  cat("\n\n")
}
## ### ATROPOS_GM_MD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_GM_MD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1067.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.55698 -0.51245  0.08047  0.54913  2.26173 
## 
## Random effects:
##  Groups   Name        Variance                       Std.Dev.          
##  subject  (Intercept) 0.0000000002327184936992856181 0.0000152551136901
##  site     (Intercept) 0.0000000000000000000000008819 0.0000000000009391
##  Residual             0.0000000000543621244491441043 0.0000073730675061
## Number of obs: 54, groups:  subject, 11; site, 3
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000757099 0.000004715   160.6
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### ATROPOS_GM_RD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_GM_RD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1067.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.69258 -0.49167  0.02313  0.54945  2.34842 
## 
## Random effects:
##  Groups   Name        Variance                Std.Dev.      
##  subject  (Intercept) 0.000000000269168850518 0.000016406366
##  site     (Intercept) 0.000000000000000007242 0.000000002691
##  Residual             0.000000000052691491108 0.000007258890
## Number of obs: 54, groups:  subject, 11; site, 3
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000695292 0.000005051   137.7
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00439263 (tol = 0.002, component 1)

p-values

p_df <- make_SL_df(models, models_0)
Plot
plot_p_SL(p_df)

Table
options(scipen = 999)
p_df %>% 
  arrange(pval_fdr)

ICC

ICC for the 35 remaining models is shown below

# extract icc_df
icc_df <- df %>% 
  filter(site != 'Hopkins') %>% 
  get_icc()

icc_df
Site
options(scipen=999)
boot_df <- readRDS(here::here('results/avg_tensor_by_roi_wide_no_Hopkins_boot_icc.rds'))
plot_ICC_site(icc_df, boot_df)

Subject
plot_ICC_subject(icc_df, boot_df)

After FDR-adjustment, 14 of the 36 variables show site effects according to the mixed models run on non-Hopkins data.

All sites (harmonized data)

Here, we report results from the mixed models carried out only on the non-Hopkins data.

models <- harmonized_df %>% 
  run_full_models()

models_0 <- harmonized_df %>% 
  run_null_models()

zero_var_model_names <- get_zero_var_model_names(models)
See models with 0 variance
for (name in zero_var_model_names){
  cat('###', name, '\n\n')
  print(summary(models[[name]]))
  cat("\n\n")
}
## ### ATROPOS_GM_AD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_GM_AD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1470.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9497 -0.6559  0.1219  0.5873  1.9042 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000018824 0.000013720
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000006596 0.000008122
## Number of obs: 74, groups:  subject, 11; site, 4
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000875530 0.000004252   205.9
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### ATROPOS_WM_AD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_WM_AD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1510.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.5691 -0.4094  0.0438  0.4148  2.2466 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000011324 0.000010641
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000003757 0.000006129
## Number of obs: 74, groups:  subject, 11; site, 4
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000891003 0.000003293   270.6
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FAST_GM_AD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_GM_AD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1610.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7396 -0.6439  0.1238  0.6573  1.9164 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000012776 0.000011303
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000005367 0.000007326
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000865351 0.000003509   246.6
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FAST_WM_AD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_WM_AD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1635.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.5691 -0.4683  0.0466  0.4440  1.9366 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000011448 0.000010700
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000003799 0.000006163
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000898376 0.000003302   272.1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FIRST_THALAMUS_AD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FIRST_THALAMUS_AD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1493.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.89528 -0.56640  0.04294  0.66129  1.63512 
## 
## Random effects:
##  Groups   Name        Variance                         Std.Dev.         
##  subject  (Intercept) 0.000000000564024297156820700447 0.000023749195716
##  site     (Intercept) 0.000000000000000000000000000121 0.000000000000011
##  Residual             0.000000000237605972500836866663 0.000015414472826
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000881333 0.000007374   119.5
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### JLF_GM_AD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_GM_AD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1581.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2746 -0.6955  0.1126  0.7008  1.6449 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.  
##  subject  (Intercept) 0.00000000043022 0.00002074
##  site     (Intercept) 0.00000000000000 0.00000000
##  Residual             0.00000000006889 0.00000830
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000918269 0.000006325   145.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### JLF_WM_AD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_WM_AD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1554.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1783 -0.5615 -0.0164  0.6371  1.9517 
## 
## Random effects:
##  Groups   Name        Variance        Std.Dev.  
##  subject  (Intercept) 0.0000000003024 0.00001739
##  site     (Intercept) 0.0000000000000 0.00000000
##  Residual             0.0000000001069 0.00001034
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000972874 0.000005375     181
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### JLF_THALAMUS_AD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_THALAMUS_AD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1499.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.91188 -0.58387  0.09154  0.69656  1.55554 
## 
## Random effects:
##  Groups   Name        Variance        Std.Dev.  
##  subject  (Intercept) 0.0000000005985 0.00002446
##  site     (Intercept) 0.0000000000000 0.00000000
##  Residual             0.0000000002165 0.00001472
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000873548 0.000007565   115.5
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### MIMOSA_LESION_AD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: MIMOSA_LESION_AD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1307.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2819 -0.6582  0.0197  0.4861  2.4428 
## 
## Random effects:
##  Groups   Name        Variance                        Std.Dev.          
##  subject  (Intercept) 0.00000001033819770882159319180 0.0001016769281048
##  site     (Intercept) 0.00000000000000000000000007073 0.0000000000002659
##  Residual             0.00000000230136384853161502135 0.0000479725322297
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept) 0.00124147 0.00003114   39.86
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### ATROPOS_GM_FA 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_GM_FA ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -615.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5514 -0.5452 -0.0296  0.4591  4.5337 
## 
## Random effects:
##  Groups   Name        Variance    Std.Dev.
##  subject  (Intercept) 0.000039622 0.006295
##  site     (Intercept) 0.000000000 0.000000
##  Residual             0.000007389 0.002718
## Number of obs: 74, groups:  subject, 11; site, 4
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.160736   0.001926   83.45
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### ATROPOS_WM_FA 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_WM_FA ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -534.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1642 -0.5512  0.0460  0.5211  2.1663 
## 
## Random effects:
##  Groups   Name        Variance   Std.Dev.
##  subject  (Intercept) 0.00043293 0.02081 
##  site     (Intercept) 0.00000000 0.00000 
##  Residual             0.00001823 0.00427 
## Number of obs: 74, groups:  subject, 11; site, 4
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.387185   0.006295   61.51
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FAST_GM_FA 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_GM_FA ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -642.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1355 -0.4601 -0.0775  0.5206  4.2706 
## 
## Random effects:
##  Groups   Name        Variance   Std.Dev.
##  subject  (Intercept) 0.00005969 0.007726
##  site     (Intercept) 0.00000000 0.000000
##  Residual             0.00001007 0.003174
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.173234   0.002357   73.48
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FAST_WM_FA 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_WM_FA ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -572.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.99637 -0.54746  0.06809  0.49398  2.41775 
## 
## Random effects:
##  Groups   Name        Variance   Std.Dev.
##  subject  (Intercept) 0.00056946 0.023863
##  site     (Intercept) 0.00000000 0.000000
##  Residual             0.00002005 0.004478
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.403270   0.007213   55.91
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FIRST_THALAMUS_FA 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FIRST_THALAMUS_FA ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -526.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.92394 -0.68760  0.09937  0.63262  2.03870 
## 
## Random effects:
##  Groups   Name        Variance   Std.Dev.
##  subject  (Intercept) 0.00036149 0.019013
##  site     (Intercept) 0.00000000 0.000000
##  Residual             0.00004206 0.006485
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.30928    0.00578   53.51
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### JLF_GM_FA 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_GM_FA ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -690.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.67793 -0.65587  0.09039  0.70079  1.91229 
## 
## Random effects:
##  Groups   Name        Variance    Std.Dev.
##  subject  (Intercept) 0.000042752 0.006539
##  site     (Intercept) 0.000000000 0.000000
##  Residual             0.000005261 0.002294
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.149733   0.001989   75.29
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### JLF_THALAMUS_FA 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_THALAMUS_FA ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -473.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.78535 -0.67347 -0.08739  0.59236  3.04956 
## 
## Random effects:
##  Groups   Name        Variance   Std.Dev.
##  subject  (Intercept) 0.00025223 0.01588 
##  site     (Intercept) 0.00000000 0.00000 
##  Residual             0.00009467 0.00973 
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.319816   0.004916   65.06
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### MIMOSA_LESION_FA 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: MIMOSA_LESION_FA ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -350.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2050 -0.3210  0.0810  0.3787  3.7290 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 0.001507 0.03882 
##  site     (Intercept) 0.000000 0.00000 
##  Residual             0.000434 0.02083 
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.32570    0.01195   27.27
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### ATROPOS_GM_MD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_GM_MD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1483.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2608 -0.5519  0.0936  0.4934  2.1900 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000022050 0.000014849
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000005204 0.000007214
## Number of obs: 74, groups:  subject, 11; site, 4
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000753768 0.000004561   165.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### ATROPOS_WM_MD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_WM_MD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1596
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9838 -0.4482  0.1488  0.5683  1.6362 
## 
## Random effects:
##  Groups   Name        Variance          Std.Dev.   
##  subject  (Intercept) 0.000000000241102 0.000015527
##  site     (Intercept) 0.000000000000000 0.000000000
##  Residual             0.000000000008663 0.000002943
## Number of obs: 74, groups:  subject, 11; site, 4
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000615973 0.000004695   131.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FAST_GM_MD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_GM_MD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1632.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3337 -0.5809  0.0119  0.4674  2.4205 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000017616 0.000013272
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000003744 0.000006118
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000736281 0.000004062   181.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FAST_WM_MD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_WM_MD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1715.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0897 -0.5025  0.1614  0.7020  1.6303 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000027296 0.000016522
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000001059 0.000003254
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000611652 0.000004995   122.4
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FIRST_THALAMUS_MD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FIRST_THALAMUS_MD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1570.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.92049 -0.48499 -0.00084  0.59701  2.01372 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000020565 0.000014341
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000008974 0.000009473
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000663841 0.000004457   148.9
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### JLF_GM_MD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_GM_MD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1605.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.00814 -0.67017  0.06203  0.68540  1.95853 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000049202 0.000022182
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000004805 0.000006932
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000802454 0.000006735   119.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### JLF_WM_MD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_WM_MD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1652.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.65999 -0.63934  0.01372  0.67020  1.81063 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000027771 0.000016665
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000002641 0.000005139
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000649447 0.000005059   128.4
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### JLF_THALAMUS_MD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_THALAMUS_MD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1569.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0016 -0.5577  0.0336  0.6803  1.8869 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000030174 0.000017371
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000008641 0.000009296
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000651171 0.000005344   121.8
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### MIMOSA_LESION_MD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: MIMOSA_LESION_MD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1314.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9024 -0.5410 -0.0192  0.4109  3.2522 
## 
## Random effects:
##  Groups   Name        Variance       Std.Dev.  
##  subject  (Intercept) 0.000000005411 0.00007356
##  site     (Intercept) 0.000000000000 0.00000000
##  Residual             0.000000002272 0.00004766
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept) 0.00091150 0.00002284   39.91
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### ATROPOS_GM_RD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_GM_RD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1484.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3266 -0.5036  0.0753  0.5078  2.2200 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000025272 0.000015897
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000005024 0.000007088
## Number of obs: 74, groups:  subject, 11; site, 4
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000692596 0.000004869   142.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### ATROPOS_WM_RD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_WM_RD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1602.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.39757 -0.52379  0.01794  0.64805  2.18877 
## 
## Random effects:
##  Groups   Name        Variance          Std.Dev.   
##  subject  (Intercept) 0.000000000407377 0.000020184
##  site     (Intercept) 0.000000000000000 0.000000000
##  Residual             0.000000000007204 0.000002684
## Number of obs: 74, groups:  subject, 11; site, 4
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000478127 0.000006094   78.45
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FAST_GM_RD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_GM_RD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1634.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5475 -0.4727 -0.0315  0.5050  2.5274 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000022075 0.000014858
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000003514 0.000005928
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000671485 0.000004531   148.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FAST_WM_RD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_WM_RD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1712.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.98135 -0.72897  0.07567  0.74622  2.14054 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000049321 0.000022208
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000001013 0.000003182
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000467961 0.000006706   69.78
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FIRST_THALAMUS_RD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FIRST_THALAMUS_RD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1604.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.71875 -0.49645 -0.00331  0.50080  2.37664 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000022539 0.000015013
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000005451 0.000007383
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000554630 0.000004605   120.5
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### JLF_GM_RD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_GM_RD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1613.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9257 -0.6679  0.0833  0.6353  2.0761 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000053911 0.000023219
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000004216 0.000006493
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept) 0.00074452 0.00000704   105.8
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### JLF_WM_RD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_WM_RD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1704.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9020 -0.6091 -0.1556  0.6481  2.1291 
## 
## Random effects:
##  Groups   Name        Variance                   Std.Dev.        
##  subject  (Intercept) 0.000000000462825961584492 0.00002151339029
##  site     (Intercept) 0.000000000000000000007618 0.00000000008728
##  Residual             0.000000000011480548228789 0.00000338829577
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000487319 0.000006498   74.99
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### JLF_THALAMUS_RD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_THALAMUS_RD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1580.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.55747 -0.62126  0.09323  0.66516  1.91966 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000028167 0.000016783
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000007436 0.000008623
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000539506 0.000005155   104.7
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### MIMOSA_LESION_RD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: MIMOSA_LESION_RD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1309.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0764 -0.5268 -0.0541  0.3975  3.4756 
## 
## Random effects:
##  Groups   Name        Variance       Std.Dev.  
##  subject  (Intercept) 0.000000005008 0.00007077
##  site     (Intercept) 0.000000000000 0.00000000
##  Residual             0.000000002472 0.00004972
## Number of obs: 80, groups:  subject, 11; site, 4
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept) 0.00074756 0.00002208   33.86
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

p-values

p_df <- make_SL_df(models, models_0)
Plot
plot_p_SL(p_df)

Table
options(scipen = 999)
p_df %>% 
  arrange(pval_fdr)

ICC

ICC for the 35 remaining models is shown below

# extract icc_df
icc_df <- harmonized_df %>% 
  get_icc()

icc_df
Site
options(scipen=999)
boot_df <- readRDS(here::here('results/avg_tensor_by_roi_wide_harmonized_boot_icc.rds'))
plot_ICC_site(icc_df, boot_df)

Subject
plot_ICC_subject(icc_df, boot_df)

After FDR-adjustment, 14 of the 36 variables show site effects according to the mixed models run on non-Hopkins data.

No Hopkins (harmonized)

Here, we report results from the mixed models carried out only on the non-Hopkins data.

models <- harmonized_no_hopkins_df %>% 
  run_full_models()

models_0 <- harmonized_no_hopkins_df %>% 
  run_null_models()

zero_var_model_names <- get_zero_var_model_names(models)
See models with 0 variance
for (name in zero_var_model_names){
  cat('###', name, '\n\n')
  print(summary(models[[name]]))
  cat("\n\n")
}
## ### ATROPOS_GM_AD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_GM_AD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1057.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.59679 -0.56995 -0.00606  0.47152  2.21326 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000019723 0.000014044
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000007148 0.000008455
## Number of obs: 54, groups:  subject, 11; site, 3
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000880281 0.000004398   200.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### ATROPOS_WM_AD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_WM_AD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1123.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.56177 -0.61044  0.04468  0.46739  1.89250 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000011658 0.000010797
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000001751 0.000004185
## Number of obs: 54, groups:  subject, 11; site, 3
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000898208 0.000003308   271.5
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FAST_GM_AD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_GM_AD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1197.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.59468 -0.54732 -0.07949  0.49328  2.18728 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000012838 0.000011331
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000005379 0.000007334
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000869450 0.000003549     245
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FAST_WM_AD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_WM_AD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1250.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.67798 -0.51885 -0.02254  0.49969  1.87665 
## 
## Random effects:
##  Groups   Name        Variance        Std.Dev.   
##  subject  (Intercept) 0.0000000001268 0.000011260
##  site     (Intercept) 0.0000000000000 0.000000000
##  Residual             0.0000000000184 0.000004289
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000904548 0.000003441   262.9
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FIRST_THALAMUS_AD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FIRST_THALAMUS_AD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1108.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.70139 -0.41705 -0.00096  0.41623  1.73233 
## 
## Random effects:
##  Groups   Name        Variance        Std.Dev.  
##  subject  (Intercept) 0.0000000007067 0.00002658
##  site     (Intercept) 0.0000000000000 0.00000000
##  Residual             0.0000000002326 0.00001525
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000881549 0.000008261   106.7
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### JLF_GM_AD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_GM_AD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1170.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0924 -0.5592 -0.1059  0.6010  1.7205 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000042429 0.000020598
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000007354 0.000008575
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000924270 0.000006312   146.4
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### JLF_WM_AD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_WM_AD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1165.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.48121 -0.50891 -0.05478  0.59726  1.96538 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000036606 0.000019133
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000008413 0.000009172
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000977541 0.000005893   165.9
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### JLF_THALAMUS_AD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_THALAMUS_AD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1111.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.74833 -0.43997 -0.02685  0.49318  1.72483 
## 
## Random effects:
##  Groups   Name        Variance        Std.Dev.  
##  subject  (Intercept) 0.0000000007461 0.00002732
##  site     (Intercept) 0.0000000000000 0.00000000
##  Residual             0.0000000002162 0.00001470
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000873240 0.000008458   103.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### MIMOSA_LESION_AD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: MIMOSA_LESION_AD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1005.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.25114 -0.66721 -0.04421  0.52226  2.38080 
## 
## Random effects:
##  Groups   Name        Variance       Std.Dev.  
##  subject  (Intercept) 0.000000012896 0.00011356
##  site     (Intercept) 0.000000000000 0.00000000
##  Residual             0.000000001061 0.00003257
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept) 0.00125981 0.00003451   36.51
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### ATROPOS_GM_FA 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_GM_FA ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -434.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9565 -0.5071 -0.0176  0.4376  3.9716 
## 
## Random effects:
##  Groups   Name        Variance    Std.Dev.
##  subject  (Intercept) 0.000047638 0.006902
##  site     (Intercept) 0.000000000 0.000000
##  Residual             0.000007877 0.002807
## Number of obs: 54, groups:  subject, 11; site, 3
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.161855   0.002118   76.41
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### ATROPOS_WM_FA 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_WM_FA ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -397.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.06094 -0.58666  0.06531  0.60149  2.30536 
## 
## Random effects:
##  Groups   Name        Variance   Std.Dev.
##  subject  (Intercept) 0.00046197 0.021493
##  site     (Intercept) 0.00000000 0.000000
##  Residual             0.00001105 0.003324
## Number of obs: 54, groups:  subject, 11; site, 3
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.390175   0.006497   60.05
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FAST_GM_FA 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_GM_FA ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -465.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8100 -0.4464 -0.0713  0.4812  3.7837 
## 
## Random effects:
##  Groups   Name        Variance   Std.Dev.
##  subject  (Intercept) 0.00006570 0.008105
##  site     (Intercept) 0.00000000 0.000000
##  Residual             0.00001136 0.003371
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.174061   0.002484   70.09
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FAST_WM_FA 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_WM_FA ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -420.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.97762 -0.50314 -0.00207  0.51873  2.26368 
## 
## Random effects:
##  Groups   Name        Variance   Std.Dev.
##  subject  (Intercept) 0.00058732 0.02423 
##  site     (Intercept) 0.00000000 0.00000 
##  Residual             0.00001832 0.00428 
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.404517   0.007329    55.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FIRST_THALAMUS_FA 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FIRST_THALAMUS_FA ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -385.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7845 -0.6547  0.1306  0.6505  2.0661 
## 
## Random effects:
##  Groups   Name        Variance   Std.Dev.
##  subject  (Intercept) 0.00041177 0.020292
##  site     (Intercept) 0.00000000 0.000000
##  Residual             0.00003988 0.006315
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.309028   0.006174   50.05
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### JLF_WM_FA 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_WM_FA ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -450.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6524 -0.6892  0.1759  0.5564  2.0740 
## 
## Random effects:
##  Groups   Name        Variance                  Std.Dev.       
##  subject  (Intercept) 0.00048879976665297715250 0.0221088164915
##  site     (Intercept) 0.00000000000000000001728 0.0000000001314
##  Residual             0.00001037093894721858281 0.0032203942223
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.428463   0.006679   64.15
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### JLF_THALAMUS_FA 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_THALAMUS_FA ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -341.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6245 -0.6431 -0.1295  0.6183  2.5472 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  subject  (Intercept) 0.0003010 0.01735 
##  site     (Intercept) 0.0000000 0.00000 
##  Residual             0.0001039 0.01020 
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 0.321027   0.005399   59.46
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### MIMOSA_LESION_FA 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: MIMOSA_LESION_FA ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -285.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1851 -0.5027  0.1120  0.4741  4.0323 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  subject  (Intercept) 0.0018111 0.04256 
##  site     (Intercept) 0.0000000 0.00000 
##  Residual             0.0002298 0.01516 
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  0.32520    0.01298   25.05
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### ATROPOS_GM_MD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_GM_MD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1061.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.74244 -0.50704 -0.00033  0.40066  2.44188 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000023279 0.000015258
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000006331 0.000007957
## Number of obs: 54, groups:  subject, 11; site, 3
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000757206 0.000004734   159.9
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### ATROPOS_WM_MD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_WM_MD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1175.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.34089 -0.42407  0.00377  0.57602  2.29103 
## 
## Random effects:
##  Groups   Name        Variance          Std.Dev.  
##  subject  (Intercept) 0.000000000242658 0.00001558
##  site     (Intercept) 0.000000000000000 0.00000000
##  Residual             0.000000000004408 0.00000210
## Number of obs: 54, groups:  subject, 11; site, 3
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000619108 0.000004706   131.6
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FAST_GM_MD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_GM_MD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1204.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.02940 -0.44962 -0.04102  0.45851  2.65693 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000017208 0.000013118
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000004435 0.000006659
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept) 0.00073933 0.00000405   182.5
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FAST_WM_MD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_WM_MD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1302.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2217 -0.6082  0.2078  0.6628  1.9767 
## 
## Random effects:
##  Groups   Name        Variance          Std.Dev.   
##  subject  (Intercept) 0.000000000273473 0.000016537
##  site     (Intercept) 0.000000000000000 0.000000000
##  Residual             0.000000000005432 0.000002331
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000615070 0.000004995   123.1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FIRST_THALAMUS_MD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FIRST_THALAMUS_MD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1163.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.70260 -0.41882 -0.09437  0.58217  1.94969 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000023837 0.000015439
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000009383 0.000009687
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000664123 0.000004825   137.6
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### JLF_GM_MD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_GM_MD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1182.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.17470 -0.61022  0.01029  0.64407  1.97465 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.  
##  subject  (Intercept) 0.00000000048376 0.00002199
##  site     (Intercept) 0.00000000000000 0.00000000
##  Residual             0.00000000005639 0.00000751
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000806894 0.000006704   120.4
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### JLF_WM_MD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_WM_MD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1227
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.35390 -0.53871  0.00628  0.50636  1.74888 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000028359 0.000016840
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000002528 0.000005028
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept) 0.00065111 0.00000512   127.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### JLF_THALAMUS_MD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_THALAMUS_MD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1155.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.80432 -0.45017  0.02758  0.58016  2.04464 
## 
## Random effects:
##  Groups   Name        Variance        Std.Dev.  
##  subject  (Intercept) 0.0000000003319 0.00001822
##  site     (Intercept) 0.0000000000000 0.00000000
##  Residual             0.0000000001045 0.00001022
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000650186 0.000005654     115
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### MIMOSA_LESION_MD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: MIMOSA_LESION_MD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1022.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.04787 -0.60690  0.02819  0.50027  2.43177 
## 
## Random effects:
##  Groups   Name        Variance        Std.Dev.  
##  subject  (Intercept) 0.0000000073893 0.00008596
##  site     (Intercept) 0.0000000000000 0.00000000
##  Residual             0.0000000008431 0.00002904
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept) 0.0009252  0.0000262   35.32
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### ATROPOS_GM_RD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_GM_RD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1059.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.74161 -0.42889  0.02642  0.47831  2.42025 
## 
## Random effects:
##  Groups   Name        Variance                       Std.Dev.          
##  subject  (Intercept) 0.0000000002688303057123714644 0.0000163960454291
##  site     (Intercept) 0.0000000000000000000000002266 0.0000000000004761
##  Residual             0.0000000000637590763813453111 0.0000079849280762
## Number of obs: 54, groups:  subject, 11; site, 3
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000695413 0.000005069   137.2
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### ATROPOS_WM_RD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: ATROPOS_WM_RD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1166
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.30250 -0.55062 -0.08627  0.62789  2.13673 
## 
## Random effects:
##  Groups   Name        Variance          Std.Dev.   
##  subject  (Intercept) 0.000000000421592 0.000020533
##  site     (Intercept) 0.000000000000000 0.000000000
##  Residual             0.000000000004868 0.000002206
## Number of obs: 54, groups:  subject, 11; site, 3
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000479304 0.000006199   77.32
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FAST_GM_RD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_GM_RD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1201.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1286 -0.4092 -0.0490  0.4789  2.7077 
## 
## Random effects:
##  Groups   Name        Variance                   Std.Dev.        
##  subject  (Intercept) 0.000000000216717387668859 0.00001472132425
##  site     (Intercept) 0.000000000000000000005284 0.00000000007269
##  Residual             0.000000000045005642021951 0.00000670862445
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000674042 0.000004525     149
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FAST_WM_RD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FAST_WM_RD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1274.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9017 -0.5845  0.0928  0.6056  2.0575 
## 
## Random effects:
##  Groups   Name        Variance          Std.Dev.   
##  subject  (Intercept) 0.000000000501142 0.000022386
##  site     (Intercept) 0.000000000000000 0.000000000
##  Residual             0.000000000008631 0.000002938
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000470076 0.000006761   69.53
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### FIRST_THALAMUS_RD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: FIRST_THALAMUS_RD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1184.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.44872 -0.64361 -0.08354  0.52868  1.98840 
## 
## Random effects:
##  Groups   Name        Variance        Std.Dev.   
##  subject  (Intercept) 0.0000000002388 0.000015452
##  site     (Intercept) 0.0000000000000 0.000000000
##  Residual             0.0000000000622 0.000007887
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000555019 0.000004772   116.3
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### JLF_GM_RD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_GM_RD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1184.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.10833 -0.54398  0.01153  0.63777  2.01162 
## 
## Random effects:
##  Groups   Name        Variance        Std.Dev.   
##  subject  (Intercept) 0.0000000005320 0.000023065
##  site     (Intercept) 0.0000000000000 0.000000000
##  Residual             0.0000000000531 0.000007287
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept) 0.00074820 0.00000702   106.6
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### JLF_WM_RD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_WM_RD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1251.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.75003 -0.46773  0.00206  0.46903  1.90570 
## 
## Random effects:
##  Groups   Name        Variance         Std.Dev.   
##  subject  (Intercept) 0.00000000046178 0.000021489
##  site     (Intercept) 0.00000000000000 0.000000000
##  Residual             0.00000000001385 0.000003721
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000487561 0.000006498   75.04
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### JLF_THALAMUS_RD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: JLF_THALAMUS_RD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1156.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3887 -0.6660  0.1540  0.6069  1.8568 
## 
## Random effects:
##  Groups   Name        Variance        Std.Dev.  
##  subject  (Intercept) 0.0000000002898 0.00001702
##  site     (Intercept) 0.0000000000000 0.00000000
##  Residual             0.0000000001041 0.00001020
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##                Estimate  Std. Error t value
## (Intercept) 0.000538321 0.000005304   101.5
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## 
## 
## 
## ### MIMOSA_LESION_RD 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: MIMOSA_LESION_RD ~ (1 | subject) + (1 | site)
##    Data: data
## 
## REML criterion at convergence: -1021.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2834 -0.7032  0.0896  0.6080  2.2461 
## 
## Random effects:
##  Groups   Name        Variance        Std.Dev.  
##  subject  (Intercept) 0.0000000069865 0.00008359
##  site     (Intercept) 0.0000000000000 0.00000000
##  Residual             0.0000000008635 0.00002939
## Number of obs: 60, groups:  subject, 11; site, 3
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept) 0.00075895 0.00002549   29.77
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

p-values

p_df <- make_SL_df(models, models_0)
Plot
plot_p_SL(p_df)

Table
options(scipen = 999)
p_df %>% 
  arrange(pval_fdr)

ICC

ICC for the 35 remaining models is shown below

# extract icc_df
icc_df <- harmonized_no_hopkins_df %>% 
  get_icc()

icc_df
Site
options(scipen=999)
boot_df <- readRDS(here::here('results/avg_tensor_by_roi_wide_no_Hopkins_harmonized_boot_icc.rds'))
plot_ICC_site(icc_df, boot_df)

Subject
plot_ICC_subject(icc_df, boot_df)

After FDR-adjustment, 14 of the 36 variables show site effects according to the mixed models run on non-Hopkins data.